##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)
library(ggpubr)
library(ggsci)
library(ggrepel)

##########################################################################################

option_list <- list(
    make_option(c("--input_file"), type = "character") ,
    make_option(c("--base_info_file"), type = "character") ,
    make_option(c("--ccf_file"), type = "character") ,
    make_option(c("--gene_list"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    input_file <- paste(work_dir,"/images/evolutionTime/timing_molecular_clock.tsv",sep="")
    ccf_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutationTime/result/All_CCF_mutTime.addShare.tsv"
    gene_list <- "~/20220915_gastric_multiple/dna_combinePublic/images/selectGCClone/GCClone_gene.all_record.list"
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/evolutionTime"

    base_info_file <- "~/20220915_gastric_multiple/dna_combinePublic/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
ccf_file <- opt$ccf_file
gene_list <- opt$gene_list
sample_info <- opt$sample_info
base_info_file <- opt$base_info_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

dat_ccf <- fread( ccf_file )
dat_input <- data.frame(fread(input_file))
dat_info <- data.frame(fread(sample_info))
dat_base_info <- data.frame(fread(base_info_file))
dat_gene <- data.frame(fread(gene_list , header = T))
colnames(dat_gene) <- "Gene_Symbol"

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

dat_ccf$Location <- paste( dat_ccf$Chr , dat_ccf$Start_Position , 
    dat_ccf$REF , dat_ccf$ALT , sep=":"  )

###########################################################################################
dat_base_info <- dat_base_info[,c("Patient" , "Alcohol" , "CNV_Type")]
dat_base_info <- subset(dat_base_info , CNV_Type!="")
dat_input <- unique(merge( dat_input , dat_base_info[,c("Patient" , "Alcohol", "CNV_Type")] , by.x = "ID" , by.y = "Patient" ))

###########################################################################################
result_driver <- c()

for( Sample in unique(dat_input$ID) ){

    tumors <- subset( dat_info , ID == Sample )$Tumor

    ## 确定driver突变
    tmp <- subset( dat_ccf , Sample %in% tumors & Variant_Classification %in% Variant_Type & Hugo_Symbol %in% dat_gene$Gene_Symbol )

    ## 判断突变的数量
    tmp_driver <- tmp %>% 
    group_by(Location) %>%
    summarize( MutTumor = length(Sample) )

    if( nrow(tmp_driver) == 0){
        class <- "NoDriver"
        tmp_res <- data.frame( Normal = Sample , DriverClass = class , 
        Share_gene = "" , Priver_gene = "" , Share_gene_variant = "" , Private_gene_variant = "") 

    }else if( nrow(tmp_driver) > 0 ){
        share_driver <- length(which(tmp_driver$MutTumor == length(tumors)))

        ## 判断是否存在部分IM与GC共享
        tmp1 <- tmp
        tmp_driver1 <- tmp1 %>% 
        group_by(Location ) %>%
        summarize( MutClass = paste0(Class , collapse = "_") )
        share_driver1 <- length(which(tmp_driver1$MutClass %in% c("IM_IGC" , "IM_DGC" , "IGC_IM" , "DGC_IM")))

        if(share_driver > 0 ){
            class <- "ShareDriver"
        }else if(share_driver1 > 0){
            class <- "PartDriver"            
        }else{
            class <- "PrivateDriver"
        }

        share_gene <- unique(subset( tmp , Location %in% subset( tmp_driver , MutTumor == length(tumors) )$Location )$Hugo_Symbol)
        private_gene <- unique(subset( tmp , Location %in% subset( tmp_driver , MutTumor < length(tumors) )$Location )$Hugo_Symbol)

        share_gene_class <- subset( tmp , Location %in% subset( tmp_driver , MutTumor == length(tumors) )$Location )$Variant_Classification
        private_gene_class <- subset( tmp , Location %in% subset( tmp_driver , MutTumor < length(tumors) )$Location )$Variant_Classification

        tmp_res <- data.frame( Normal = Sample , DriverClass = class , 
        Share_gene = paste0(share_gene , collapse = ",") , Priver_gene = paste0(private_gene , collapse = ",") ,
        Share_gene_variant = paste0(share_gene_class , collapse = ",") , Private_gene_variant = paste0(private_gene_class , collapse = ",") 
        )
    }

   
    result_driver <- rbind( result_driver , tmp_res )
}

###########################################################################################

dat_input2 <- merge( dat_input , result_driver , by.x = "ID" , by.y = "Normal" )
dat_input2$plotType <- ifelse( dat_input2$Share_gene != "" , "TrunkDriver" , "Other" )
dat_input2$plotType <- factor( dat_input2$plotType , levels = c("TrunkDriver" , "Other") , order = T )
dat_input2$Class <- factor( dat_input2$Class , levels = c("IGC" , "DGC") , order = T )

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

col_tmp <- c(
    rgb(red=179,green=60,blue=59,alpha=255,max=255) ,
    rgb(red=14,green=90,blue=170,alpha=255,max=255)
)

###########################################################################################
## 比较GS和CIN
for( class in c("All" , "IGC" , "DGC" , "TrunkDriver" , "Other" ,  "Drink" , "No" ) ){

    if(class == "All"){
        dat_input3 <- dat_input2
    }else if(class %in% c("IGC" , "DGC")){
        dat_input3 <- subset( dat_input2 , Class == class )
    }else if(class %in% c("TrunkDriver" , "Other")){
        dat_input3 <- subset( dat_input2 , plotType == class )
    }else if(class %in% c("Drink" , "No")){
        dat_input3 <- subset( dat_input2 , Alcohol == class )
    }

    p <- wilcox.test(subset(dat_input3 , CNV_Type=="GS")$mean , subset(dat_input3 , CNV_Type=="CIN")$mean)$p.value

    if( p < 0.01 ){
        p_text <- trans(p)
    }else{
        p_text <- paste0( "P == " , round(as.numeric(p) , 3) ) 
    }

    dat_input3$p_text <- ""
    dat_input3$p_text[1] <- p_text
    dat_input3$label <- ifelse( dat_input3$Share_gene != "" , paste0( dat_input3$Share_gene , "," , dat_input3$Priver_gene ) , "" )

    y_lab <- "Number of years for \nIM to GC progression"
    y_max <- 20 + 1

    plot <- ggplot( dat_input3 , aes( x = CNV_Type , y = mean , color = CNV_Type ) ) +
        geom_boxplot(size = 1.2 , outlier.shape = NA ) + ## 去除散点，加粗线
        scale_color_manual(values=col_tmp) +
        scale_fill_manual(values =col_tmp) +
        geom_jitter(position=position_jitter(0.2 , seed = 1),aes(color=CNV_Type)) +
        geom_text(aes(label=p_text , y = y_max , x = 1.5),parse = TRUE,size=5 , color = "black", face='bold') +
        xlab(NULL) +
        ylab(y_lab)+
        theme_bw() +
        # geom_text_repel( position = position_jitter(seed = 1) ,
        #        aes(label = label) , size = 2 , col = 'black' , face = 'bold') +
        theme(
            legend.position = 'none',
            legend.title = element_blank() ,
            panel.grid.major=element_blank(),
            panel.grid.minor=element_blank(),
            panel.background = element_blank(),
            panel.border = element_blank(),
            plot.title = element_text(size = 12,color="black",face='bold'),
            legend.text = element_text(size = 12,color="black",face='bold'),
            axis.text.y = element_text(size = 12,color="black",face='bold'),
            axis.title.x = element_text(size = 12,color="black",face='bold'),
            axis.title.y = element_text(size = 14,color="black",face='bold'),
            #axis.text.x = element_text(size = 14,color="black",face='bold',angle = 45,hjust = 1) ,
            axis.text.x = element_text(size = 14,color="black",face='bold') ,
            axis.ticks.length = unit(0.2, "cm") ,
            strip.text.x = element_text(size = 17, colour = "black",face='bold') ,
            axis.line = element_line(size = 0.5)) 

    width <- 3/1
    height <- 4.47/1
    out_name <- paste0(out_path , "/CompareGS_CIN.Time." , class , ".pdf")
    ggsave(file=out_name,plot=plot,width=width,height=height)
}

out_name <- paste0(out_path , "/CompareTrunkSMG.Time.tsv")
write.table(dat_input2 , out_name , sep = "\t" , quote = F , row.names = F )